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求解 大 型 矩阵 特征 值 问题 的 并 行 精 化 Davidson 方法 
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摘 E 针对 共享 主 存 的 并 行 计算 环境 和 微机 网 络 并 行 计 算 环 境 ， 本 文 给 出 了 求解 大 型 稀 玻 对 称 算 阵 的 部 
分 极端 特征 对 的 并 行 精 化 Davidson 方 法 ， 分 析 了 该 法 的 内 在 并 行 性 。 各 处 理 器 利用 挎 阵 的 行 块 
和 投影 子 空间 的 正 交 基 所 组 成 托 阵 的 行 块 进行 运算 ， 结 合 重新 启动 策略 求解 矩阵 多 个 特征 对 的 近 
似 值 ， 并 用 以 计算 某 型 号 机 哟 的 固有 频率 ， 在 微机 网 络 并 行 计 算 环 境 和 拥有 共享 主 存 并 行 计算 环 

， 境 IBM-P650 上 进行 了 数值 试验 。 
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1 引言 


在 科学 和 工程 技术 领域 中 ， 经 常 需要 计算 大 型 稀 踊 对 称 和 矩阵 的 若干 个 特征 对 。 解 决 这 一 问 
题 的 有 效 方法 是 正 交 投影 方法 ， 包 括 子 空间 迭代 法 ，Lanczos 方法 ，Arnoldi 方法 ，Davidson 方 
法 和 Jacobi-Davidson 方法 中 3。 为 进一步 提高 Lanczos 方法 和 Davidson 方 法 的 可 靠 性 ，Under- 
wood 提 出 了 块 Lanczos 方 法 回 ， 戴 等 提出 了 预 处 理 块 Lanczos 方 法 由 ，Crouzeix 等 人 提出 了 
块 Davidson 方法 [5.61 。 

用 正 交 投影 方法 计算 特征 值 问题 

AT = Àz (1) 
的 部 分 特征 对 时 ，Ritz 向 量 2; 可 能 收敛 很 慢 其 至 发 散 ， 针 对 这 一 问题 ， 费 提出 了 保 
留 Ritz 值 A;， 用 满足 


IA- ADs = min _ IA- Ad)alla (2) 


的 精 化 向 量 2; 代替 Ritz 向 量 g AR), Hee TT 
I(A — AsZ)Zille = omin(AVi — NW), (3) 
SO ŝi = Vidi, 2i RI (A 一 Xi 丰硕 的 最 小 奇异 值 cmin 对 应 的 右 奇 异 向 量 ， 晶 旧名 2 = 1。 即 


VI(A— NAD)" (A — AD)Vi8i = OF in Zi (4) 
个 难 证 明 $; = Vids EXIT: Ritz 值 的 使 本 次 迭代 的 残 量 范 数 |(4 —AD Fillo 最 小 的 向 量 ， 此 时 
的 残 量 范 数 为 cmin， 特 征 值 的 近似 可 取 为 Ai = 2E Aa; 
大 型 矩阵 特征 值 问题 的 计算 是 一 项 很 耗 时 的 工作 ， 尽 管 人 们 采用 多 种 方法 降低 计算 量 ， 但 
终 因 串 行 思想 的 限制 难以 有 更 大 的 提高 ， 采 用 并 行 计 算是 一 条 必由之路 。 多 人 年来， 很 多 科技 工 
Wc RA: 2007-11-13. ESA: FIRA (1962 年 1 月 生 )， 男 ， 博 士 ， 副 教授 . 研究 方向 ， 数 值 代数 . 
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作者 对 特征 值 问题 的 并 行 计算 进行 了 大 量 的 研究 。Balle 和 Cullum、Nool 和 Van der Ploeg 4} 
别 讨论 了 并 行 Lanczos 方 法 和 并 行 Jacobi-Davidson 方法 0314， 取 得 了 满意 的 效果 。 为 提高 求 
解 大 型 矩阵 特征 值 问题 块 Davidson 方法 的 计算 效率 ， 本 文 提出 并 行 精 化 块 Davidson 方 法 ， 并 
分 别 在 微机 网 络 并 行 环境 和 并 行 计算 机 IBM-P650 上 进行 了 数值 试验 。 


2 精 化 块 Davidson 方法 


WA A n BYERS BRIE, Davidson 方法 利用 Galerkin-Ritz 投影 技术 ， 计 算 矩 阵 的 Ritz 值 
和 Ritz 向 量 ， 使 用 预 处 理 技术 拓 广 投影 子 空间 V 的 基 媚 ， 在 新 的 迭代 步 中 ，Davidson 方 法 
充分 利用 已 经 获得 的 信息 扩充 新 的 子 空间 ， 因 而 在 很 少 的 迭代 步 内 得 到 满足 精度 要 求 的 
解 。Davidson 方法 通常 和 重新 开始 技术 结合 使 用 。 

算法 1 (H Davidson 方法 回 ) 

1) ”选取 正 整 数 m 和 1， 列 正 交 规范 矩阵 态 = [vi, v2,--- ,v1]， 精 度 e; 

2) 对 大 = 1,2,:…， 执 行 以 下 各 步 ; 3) Wk = AVe: 4) Hp = Vi We; 

5) 计算 Hs 的 1 个 最 大 (或 最 小 ) 特 征 对 Oni, Yki) i = 1,… ls 

6) WH Ritz a. = VkYk i t= 1, ,1l; 

7) IRIRE Tki = nal — Aiki i = 1 sh 

8) ”收敛 性 检查 ， 若 ||rgpijl2 <e, i = 1,… ,1!， 则 停止 ， 否 则 转 9); 

9) ”求解 校正 方程 Ckitki = Tei, T= 1, ,1; : 

10) #dim(V) +l<m, WW Vier = MGS(Viyth,1tk,2,… tes)» 723), AMF 11); 

11) Vi = MGS(%e1,-°+ Ee, tk, tkt tet)» 转 2) 重新 开始 。 

算法 1 中 的 MGS 代表 修正 的 Gram-Schmidt 正 交 化 过 程 ， 其 数值 稳定 性 优 于 经 典 的 Gram- 
Schmidt 正 交 化 过 程 ， 假 设 Ve = (v, ,vs) 列 正 交 规范 ，w; = tk ys, J=St+l,---,8+1, M 
算法 1 中 的 了 + = MGS(W， 大 ,1 大 ,2 大) 的 过 程 如 下 。 


fori=1toldo 
for j = s + 1 to s +1 do 
Uj = Uj (vj, Vi) Ui 
endfor 
endfor 
for i = s + 1 to s +l do 
vi = vi/|lvill2 
for j =i + 1 to s +L do 
Uj = UZ — (vj, Vi) Vi 
endfor 
endfor 


SB k UR RAR AY T ah IE E RE Chi kal — 4 的 近似 ， 4 Opi EU Akil 一 4 时 ， 校 正方 
程 Cpstki = tee (i = 1,… DNR BE, (i = 1,… ,人 在 这 种 情况 下 将 tri (i = 
1,… ,中 加 入 到 基 中 会 使 新 的 基底 线性 相关 ， 为 克服 这 个 缺点 ， 可 使 用 Olsen 预 处 理 45， 此 时 
校正 方程 为 


Chk itki = Tks + Eiki te LTes, i=1,...,l, (5) 
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由 万 ;2 = 0 得 到 
Ei = -21 OF rei/ LRC Ek, i=1,.…,l, (6) 


这 样 经 过 MGS 过 程 后 得 到 的 子 空间 的 基 标 准 正 交 。 
将 精 化 思想 应 用 于 块 Davidson 方法， 就 和 x,i 而 言 ， 


VZ (A 一 de al)7 (A 一 deil)Vi = Wa We 一 2Xk Hp 十 Nad, 


假设 


Ski = WEWe — 2iki Hk + AGI, 1 (7) 


得 到 精 化 块 Davidson 方法 ， 以 算法 2 表示 。 
算法 2 (HALE Davidson 方法 ) 
1) 选取 正 整数 mm 和 1， 列 正 交 规范 矩阵 全 = [v1,v2,… ul, Mee: 
2) 对 == 1,2,… ,执行 以 下 各 步 ; 3) Wk = AV; 4) Hx = Vi Wr; 
5) 计算 Hr 的 ! 个 最 大 (或 最 小 ) 特征 对 Oki, yea), i= 1,… ,1 
6) 计算 Ritz ARG. = Veyei, t= 1,--- hs 
7) EARE rki = (Akil — A)Žki T= 1 ls 
8) KAHERE, F irkall < s i=1, ,4， 则 停止 ， 否 则 转 9); 
9) 4; dim(V;,) 十 1<m， 求 解 Ck itki = Tki 1 = 1,---,l, Vk} = MGS(T tk,1, tee sth) 
转 3); 
10) tr dim(Vk) 十 1 > m 计 算 Ski = WE Wk 一 We Hey 十 X27, 7 = 1,.…: ,的 最 小 特征 
对 (o? ;， ĉi), 1 一 1 
11) WEAR a. = Ve, =1 1 
12) Vi = [êk1, feu], 2). 


3 ”并 行 精 化 块 Davidson 方法 


3.1 精 化 块 Davidson 方法 的 并 行 计 算 

将 ne 阶 方 阵 和 4 按 行 划分 ， 使 得 各 处 理 机 平均 分 配 和 A 的 行 数 ， 若 nw 不 是 处 理 机 数 p 的 整数 
沼 ， 则 前 间 的 处 理 机 比 后 面 的 多 拥有 一 行 。 假 设 处 理 机 Pi 中 存放 矩阵 的 部 分 为 4;，Ai 共 
Emir, BRÈ m TEM ADB i TBAB i + my 一 1 行 ， 根 据 算法 2 得 到 并 行 精 化 
块 Davidson 算法 。 

算法 3 (并 行 精 化 块 Davidson 方法 ) 

1) AREH mH PEMER V = [v1,v2,… ui), Be: 

2) Xk=1,2,---, WTU FEF; 3) Pitt Wei = AiVh; 

4) Pw Hri = Vg Wpais 5) PiP Hk = 1 Hkg: 

6) Pitt Hr 的 ! 个 最 大 (或 最 小 ) 特征 对 (Akg Yk) j= lsc ols 

7) 并 行 计算 Ritz 向量 Ek,j = Vive, j=l hs 

8) FPF BRE rey = (Arj — A)zpj, j= loak 

9) KARE, E llrejle < es 了 = 1,… MIE, AUH 10); 

10) #4 dim(V;) +l<m, 解 Cr jth =Tk j, j=l, b Veo = MGS (Vk, tk; stk)» 

tE 3); 


11) # dim(V,) +1 > m 则 PP 计算 矩阵 WiEWki， #32) WE Ww, AT 


Sk,j = We We — 2x He + X 5, j=1,-:-,l 


? 


计算 Ska = 1,… ,0) 的 最 小 特征 对 (j2) j = 4，… ds 

12) 并 行 计算 精 化 向 量 人 各 = 你 二 j=l, ,ls 

13) Vi = [êk Fea], #2). 

3.2 并行 特性 分 析 及 数据 分 配 

算法 3 的 计算 过 程 主要 是 矩阵 和 向 量 乘 ， 因 此 具有 和 良好 的 并 行 性 ， 其 中 6) 计 算 低 阶 矩 
M Ay 的 特征 对 及 11) 中 最 小 奇异 对 的 运算 量 较 小 ， 各 处 理 机 可 分 别 从 事 该 过 程 的 计算 。 

We 一 AV; 及 A, = VW: 的 并 行 计 算 。 P; 中 Wri = AiVk， 从 而 P; 可 以 独立 计算 Wk 的 Mi 
{Te BVE = (VE VET PA PRR ARTAR, EKARA AE 
EVE 的 列子 块 和 Wi 的 行 子 块 的 乘积 ， 处 理 机 PB 计算 Hri = ViTWki， 整 体 求 和 或 主 节点 机 
求 和 得 到 


P 
Hr = 2 Hrs» 
j=1 


各 节点 机 (或 主 节点 机 ) 计算 Ex 的 全 部 特征 对 ， 确 定 He 的 ! 个 特征 值 Xkj (7 = 1,… 0, th 
主 节 点 机 计算 ， 必 须 将 结果 发 送 给 各 子 进程 。 
WR Wi 的 并 行 计算 。P; AWE, Wri AR 


Sk j = WE We — gj He + Xg jI, j=l1,:--,l, 


主 节 点 机 或 各 节点 机 同时 计算 Srg 的 最 小 特征 对 (02,5, 2)), j= 1,… ,1。 
Ritz 向 Ze = 1 D RATH. Bit 


Yk 二 (Yk: Yki) Xk = (Ëk ) £1); 
HEX, 行 分 块 ， 行 分 块 的 情况 同和 的 行 分 类， 那么 Xe = (XE XE), W 
Xki = Vk, i Yk, i= 1,- … yD, 


{HEU ADSWL P; 独立 计算 Xk,i， 同 样 的 方法 用 于 并 行 计算 精 化 向 量 人 Zk,j, 7 = 1 Lo 
reg = ApjEpj 一 WeyegG = 1,… DMFT HR. BBR = (re Tet)» 
同 Xe 一 样 ， 将 Ri 行 分 块 ， 即 Re = (RE Rip) > BF 


Rk = Xkdiag(Xk 1， trey dx.) 一 WzYk, 


JARDA Ree = Xidiag (N41 Ant) —WriYer Rei 的 计算 是 独立 的 。 
Cjtkji = Thi(j = 1,… D HERI. Crj = reg 一 MA)， 若 取 MA 为 4 的 对 角 阵 
即 MA = Da, Pi” trg 的 计算 过 程 为 


(tk,j)s = (rr j)s/ Ùr — ass) s=, s itm, j=l,---,b, 


B Tk = (tket) UWP, PAA Tr 的 一 个 行 数 为 mi 的 行 块 。 
经 过 上 述 各 步 ， 各 结 点 机 中 可 得 到 Vig 的 相应 行 块 ， 然 后 进行 相应 的 正 交 化 过 程 。 
3.3 正 交 化 过 程 的 并 行 计算 
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算法 3 中 1 各 列 的 正 交 化 过 程 ， 采 用 行 分 块 的 形式 ，P; 负责 Vier 的 ms 行 的 运 
6. Veo 各 列 向 量 的 范 数 以 及 各 向 量 和 其 后 面向 量 内 积 的 计算 ， 则 由 各 处 理 机 分 别 负 
担 一 部 分 ， 为 减少 通讯 ， 将 经 典 的 Gram-Schmidt (CGS) 正 交 化 过 程 和 MGS 过 程 相 结合 ， 
记 为 CMGS 过 程 。 CMGS 过 程 的 一 部 分 环节 使 用 CGS 过 程 ， 男 一 部 份 使 用 MGS 过 程 。 
当 m -dim(W%)>1 时 ， 使 用 并 行 CGS 方 法 ， 对 处 理 机 中 的 tk,y(j = 1,… ,) 相 应 的 部 分 进 
行 正 交 修 正 ， 该 过 程 只 需 一 次 整体 求 和 ， 然 后 对 修正 过 的 ty(j = 1,… ,站 白 身 进行 MGS 标准 
正 交 化 过 程 。 


4 ”数值 试验 及 结论 


假设 p 代 表 节 点 机 的 数量 ， 如 代表 使 用 p 个 节点 机 计算 的 时 间 (单位 : BD), Sp = ti/tp R 
示 加 速 比 ，7 表 并 行 效率 (以 百分比 表示 )， 和 迭代 精度 s = 10-5。 数 值 试 验 在 IBM-P650 利 
微机 网 络 并 行 计算 环境 上 进行 ， 微 机 网 络 并 行 环境 由 安装 MPI 的 PH300 微 机 组 成 ， 内 
存 64M。IBM-P650 是 拥有 8 个 1.45GHz 的 Power4CPU、16GB 内 存 、FastT 磁 盘 阵 列 的 并 行 
计算 机 ， 在 IBM-P650 并行 机 上 用 Fortran 编译 器 编制 使 用 共享 内 存 进 行 通信 的 MPI 程 序 。 

算 例 1 nn 阶 和 矩阵 4 = (aij) 的 元 素 定 义 为 i = 了 时 ai = i ji -—j[xwif Hi<j Na, = 
Alil, WEER a= 0。w 为 4 的 半 带 宽 。 本 算 例 取 n = 7000, w= 262，A = 0.758, # 
两 种 并 行 环境 下 ， 使 用 精 化 块 Davidson 方法 和 非 精 化 块 Davidson 方法 并 行 计算 4 的 5 个 最 小 
特征 对 ， 结 果 见 表 1-3。 


Rl: 算 例 1 的 前 5 个 最 小 特征 值 


一 4.0931326 x 107? 
0.5804710 
1.0164097 
1.7284262 
2.2801648 


a FF Ww NH Fle 


表 2: 微机 并 行 环境 下 并 行 计算 算 例 1 的 情况 


并 行 精 化 块 Davidson 方法 并 行 块 Davidson 方法 


t Sp n p t Sp 
22.45 1 100 1 304.13 1 
14.48 1.55 77.5 2 187.1 1.63 


表 3: IBM-P650 上 并 行 计算 算 例 工 的 情况 
并 行 精 化 块 Davidson 方法 并 行 块 Davidson 方法 
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算 例 2 ”使 用 并 行 精 化 块 Davidson 方 法 和 非 精 化 块 Davidson 方 法 ， 并 行 计算 如 图 1 所 示 机 
REAMER., R K A AM 代表 结构 的 刚度 敌阵 和 质量 矩阵 ， 和 是 特征 值 问题 天 rz = 
AMz 的 特征 值 ， 则 结构 相应 的 固有 频率 f = VA/2r， 将 广义 特征 值 问 题 转变 成 By = iy, E 
"H B =L ML-T, K=LLT, X=1/d, RENI tki = Coles PORER Ons — B)? ME 
似 


Crt = (1+ 5-8). 
i 


和 ki 
在 IBM-P650 上 的 计算 情况 如 表 4 所 示 。 


图 1: 某 型 号 机 标的 上 翼 面 


4&4: IBM-P650 上 并 行 计算 算 例 2 的 情况 


并 行 精 化 块 Davidson 方法 并 行 块 Davidson 方法 


计算 结果 表明 ， 算 例 1 中 非 精 化 方法 和 精 化 方法 的 选 代 次 数 分 别 为 101 次 和 6 次 ， 精 化 方法 
比 非 精 化 方法 迭代 次 数 减 少 了 许多 ， 节 省 了 计算 时 间 ; 算 例 2 的 预 处 理 矩 阵 的 修正 效果 非常 
好 ， 两 种 方法 的 迭代 次 数 都 是 5 次 ， 精 化 方法 的 运行 时 间 比 非 精 化 的 稍 多 ， 加 速 比 稍 低 ， 这 是 
由 于 在 本 算 例 中 两 种 方法 都 重新 启动 一 次 ， 精 化 方法 需要 额外 计算 Skj (7 = 1,… 1) 及 其 相应 
的 最 小 特征 对 (of 5,2) (7 = 1,… ,人 增加 了 计算 量 和 相应 的 通讯 开销 。 精 化 方法 的 迭代 次 数 
不 会 超过 非 精 化 方法 的 迭代 次 数 ， 其 减少 计算 时 间 的 原因 在 于 其 可 能 减少 迭代 步 。 由 于 涉及 到 
的 运算 基本 上 都 是 矩阵 向 量 乘 ， 并 行 计算 的 效果 是 明显 的 。 
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Abstract: A parallel refined Davidson method is presented for large sparse symmetric eigenvalue 
problems, its implementation on PC network parallel systems and shared memory computer systems, 


as well as the parallelism are analyzed. The rows of the matrix are distributed to each processor, and 


the individual processors run the program by using the orthogonal basis of projection subspace and 


the rows of the matrix. Combine with the restarting scheme, the method can solve several eigenvalues 


and its associated eigenvectors of the matrix, and the method is successfully used for computing the 


frequency of a plane wing. The numerical experiments are done on the PC network parallel system 
and shared memory parallel system IBM-P650. 


Keywords: parallel computing; eigenvalue problem; Davidson method; refined method 


